#include <cuda.h>
#include <stdio.h>



#define OLAMDA (0)


extern "C"  void gettv_gpu(float *x_recon,float *x0, int width, int height, float lamda,int iterNum);
cudaArray* initCudaArrayTexture4Grad(texture<float, 3,cudaReadModeElementType> &texArray,int width,int height, int length, float *HostData);
__global__ void kernelGradDesSlice(float * Dx,int size,float mu,int indexZ);



  //texture<float, 3,cudaReadModeElementType> texPriVolume;
       texture<float, 3,cudaReadModeElementType> texVolumeART;
       texture<float, 3,cudaReadModeElementType> texVolumeGrad;
	   
	   
	   

void checkCUDAError(const char *msg)
{

    cudaError_t err = cudaGetLastError();
    if( cudaSuccess != err) 
    {
        fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString( err) );
        exit(EXIT_FAILURE);
    }                         
}

void UpdateVolume(float *tempDx, float *x_recon, int size, float t)
{

  for(int index =0;index<size;index++)
  {
    x_recon[index] = t*tempDx[index] + x_recon[index]; 
  }
}


__global__ void kernalGetObjFun(float *tmpResutls,int size,int indexZ)
{


	
	  int z= (indexZ);
	int x = (threadIdx.x) + (blockIdx.x)* 16;
	int y = (threadIdx.y) + (blockIdx.y)* 16; 

    float sumTV = 0;


    float GH = 0;
    float GV = 0;
	//float GHDIF = 0;
	//float GVDIF = 0;

    if(x> (size -1) || y> (size -1)) return;
	if( ( y<(size -1) ) && ( x<(size -1) ) )
	{
		//|x|TV
		
		GH = - tex3D(texVolumeGrad,x,y,z) + tex3D(texVolumeGrad,x  ,y+1,z);
		GV = - tex3D(texVolumeGrad,x,y,z) + tex3D(texVolumeGrad,x+1,y  ,z);


		sumTV += sqrtf(GH*GH + GV*GV);
	
		
	}

	tmpResutls[x + y*size] = sumTV;
     
}




extern "C"  void gettv_gpu(float *x_recon,float *x0, int width, int height, float lamda,int iterNum)
{

       int numVoxels = width * height;
	   int length =1;
	   
	   memcpy((void *)x_recon,(void*) x0,sizeof(float)*width*height);
	   
        float *DtmpResults = NULL;
	 
	 cudaMalloc((void **) &(DtmpResults), numVoxels*sizeof(float)); checkCUDAError("malloc objFun tmpResults");
	   
	
	   
	
	   
	  
       cudaArray *cudaVolumeGrad = initCudaArrayTexture4Grad(texVolumeGrad,width,height,length,x0);
	   
	  
	    


	   
	   //dx
       float *Dx;
	  cudaMalloc((void **) &(Dx), numVoxels*sizeof(float)); checkCUDAError("malloc gradient");
	   cudaMemset(Dx, 0, numVoxels*sizeof(float));  checkCUDAError("cudaMemset ");
	     
		 
      
	    //
	   dim3 grid_dbp(width/16, height/16);
	   dim3 testgrid(16, 16);
	   kernalGetObjFun<<<grid_dbp,testgrid>>>(DtmpResults,width,0);
	   
	  
	   
	   cudaMemcpy(x_recon,DtmpResults,sizeof(float)*numVoxels, cudaMemcpyDeviceToHost);
	   

	  
	  cudaUnbindTexture(texVolumeGrad);


	 
	  cudaFreeArray(cudaVolumeGrad);
	  
	  cudaFree(Dx);
	cudaFree(DtmpResults);

  
}



 cudaArray* initCudaArrayTexture4Grad(texture<float, 3,cudaReadModeElementType> &texArray,int width,int height, int length, float *HostData)
{
       
       cudaArray *cudaV;
	   cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();  
       cudaExtent extent;     
       extent.width= width;     
       extent.height=height;      
       extent.depth=length;   
       
	   //alloc cuda array
       cudaMalloc3DArray(&cudaV,&channelDesc,extent);
	   checkCUDAError("cudaMalloc3DArray cudaPriVolume");

	   //bind to texture and init
	   texArray.addressMode[0] = cudaAddressModeClamp;
	   texArray.addressMode[1] = cudaAddressModeClamp;
	   texArray.addressMode[2] = cudaAddressModeClamp;
	   texArray.filterMode = cudaFilterModePoint;
	   texArray.normalized = 0; 
       cudaBindTextureToArray(texArray, cudaV,channelDesc); 
	   checkCUDAError("cudaBindTextureToArray texVolume");
	   
	   cudaMemcpy3DParms copyParams = {0};
	   copyParams.srcPtr = make_cudaPitchedPtr((void*)(HostData), 
						width*sizeof(float), 
						height, 
						width);
	   copyParams.dstArray = cudaV;
	   copyParams.extent = extent;
	   copyParams.kind = cudaMemcpyHostToDevice;

	   cudaMemcpy3D(&copyParams);
	   checkCUDAError("cudaMemcpy3D cudaDynaIterCopyVolumeHostToDeviceArray");

	   return cudaV;

}